UNCLASSIFIED 


AD  NUMBER 


AD803195 


NEW  LIMITATION  CHANGE 
TO 

Approved  for  public  release,  distribution 
unlimited 


FROM 

Distribution  authorized  to  U.S.  Gov't, 
agencies  and  their  contractors;  Critical 
Technology;  NOV  1966.  Other  requests  shall 
be  referred  to  Air  Force  Weapons  Lab . , 
AFSC,  Kirtland  AFB,  NM. 


AUTHORITY 


AFWL  ltr,  30  Nov  1971 


THIS  PAGE  IS  UNCLASSIFIED 


AFWL-TR-66-141,  Vol.  I 


AFWL-TR 
66-141 
Vol.  I 


us 

OS 

r**i 

CO 

© 


« 1 


THEORETICAL  CALCULATIONS  OF  THE 
PHENOMENOLOGY  OF  HE  DETONATIONS 

Volume  I 


William  A.  Whitalcor,  Captain,  US'S? 
Edmund  A.  Nawrockt,  Captain,  USAF 

Charles  F.  Nctdhar/. 

William  W.  Troutman 


TECHNICAL  REPORT  NO.  AFWL-TR-66-141,  Vol.  I 
Nevombar  1966 


AIR  FORCE  WEAPONS  LABORATORY 
Research  and  Technology  Division 
Air  Fore*  Systems  Command 
Kirtiand  Air  Fere*  But* 

New  Mexico 


*  9 


803195 


»s«sj .»-» 


AFWL-TR-66-141,  Vol .  I 


THEORETICAL  CALCULATIONS  OF  THE  PHENOMENOLOGY  OF  HE  DETONATIONS 

Volume  I 


William  A.  Whitaker,  Captain,  USAF 
Edmund  A.  Nawrocki,  Captain,  USAF 
Charles  E.  Needham 
William  W.  Troutman 


TECHNICAL  REPORT  NO.  AFWL-TR-66-141,  Vol.  I 


This  document  is  subject 
to  special  export  controls  and 
each  transmittal  to  foreign 
governments  or  foreign 
nationals  may  be  made  only 
with  prior  approval  of  AFWL 
(WLRTH) ,  Kirtland  AFB,  N.M, 
Distribution  of  this  document 
is  limited  because  of  the 
technology  discussed. 


AFWL-TR-66-1U,  Vol,  I 


FOREWORD 


This  research  was  performed  under  Program  Element  7. 60. 08. 01. D,  Project 
5710,  Subtask  1.027,  and  was  funded  by  the  Defense  Atomic  Support  Agency  (DASA) . 

Inclusive  dates  of  research  vere  1  September  1965  to  1  June  1966.  The 
retort  was  submitted  27  October  1966  by  the  Project  Officer,  Captain  William  A. 
Whitaker,  (WLRTH) . 


f 


t 


This  report  has  been  reviewed  and  is  approved. 


WILLIAM  A.  WHITAKER 
Captain,  USAF 
Project  Officer 


RALPH  H.  PENNINGTON  X 

Colonel,  USAF 

Chief,  Theoretical  Branch 


Colonel,  USAF 

Chief,  Research  Division 


1 


1 


ii 


AFWT -TR-66-! Ai ,  Vol.  I 


ABSTRACT 


The  phenomenology  of  two  atmospheric  high-explosive  detonations  were  calculated 
theoretically.  The  first  was  a  20-short-ton  spherical  charge  of  TNT  (loading 
density — 1.56  gms/cc) .  The  second  was  a  methane-oxygen  mixture  (mole  ratio  1  to 
1.5)  contained  In  a  55-ft-radlus  balloon.  Both  detonations  took  place  at  an 
altitude  of  670  meters  (ambient  pressure  13.6  pal)  with  a  reflecting  surface 
85  feet  below  burst  point.  The  calculations,  taken  out  to  300  milliseconds  after 
detonations,  were  performed  by  using  SAP,  a  one-dimensional  Lagrangian  hydro- 
dynamic  code  and  SHELL-OIL,  a  two-dimensional  pure  Eulerlan  hydrodynamic  code. 
Volume  II  of  this  report  contains  the  details  of  the  results  In  graphical  form. 
Included  are  pressure  and  density  contours,  velocity  vector  plots,  and  wave  forms 
for  19  test  stations. 
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*  SECTION  I 

INTRODUCTION 

The  Air  Force  Weapon*  Laboratory  la  engaged  In  a  aerie*  of  theoretical 
calculation*  of  the  phenomenology  of  atmoepheric  nuclear  weapon  detonation*. 

These  calculation*  are  based  chiefly  upon  a  series  of  large  compute*  codes,  to 
enable  prediction  from  essentially  first  principles,  of  the  phenomaoology  from 
microseconds  to  minutes  after  a  detonation.  The  purpose  of  this  undertaking  is 
to  provide  a  best  theoretical  pict*"-*  of  the  phenonenoiogy  of  atmospheric  deto¬ 
nations  for  use  In  radar  and  re-entry  vehicle  vulnerability  studies. 

Three  main  codes  are  used  In  this  program:  SF UTTER,  a  one-dlmenslonal 
radiation  transport  Lag r mg lan  hydrodynamic  code;  SAP,  a  one-dlmenslonal 

*  Lag rang fan  hydrodynamic  code;  and  SHELL,  a  two-dimensional  Euler lan  hydrodynamic 
code.  The  first,  SPUTTER,  takes  the  radiative  output  directly  from  weapons 

*  design  calculations,  deposits  the  energy  in  air,  and  calculates  the  redlative 
and  hydrodynamic  growth  of  the  fireball  to  a  time  of  about  one  second.  It  uses 
a  multi-frequency  transport  scheme  involving  typically  20  frequency  groups.  The 
output  of  this  code  is  used  aa  Input  to  SAP,  which  calculates  the  hydrodynamic 
expansion  of  the  shock  wave  ct  various  angles  to  the  horizontal .  This  calculation 
can  proceed  to  great  distance*  so  that  even  the  shock  on  the  ground  can  be 
estimated,  but  It  is  chiefly  importt  t  In  the  strong-shock  region  of  interest 

In  R/V  vulnerability.  The  output  of  SPUTTER  la  also  used  as  input  to  SHELL, 
which  calculates  fireball  growth  and  rise  to  late  times.  The  output  of  chess 
codes  Is  processed  to  give  the  desired  observables  stub  as  pressure,  temperature, 
and  velocity. 

The  first  shot  so  calculated  u  BLUECILL,  a  high-altitude  detonation  of  the 
DOMINIC-FISHBOVL  nuclear  test  series.  This  particular  shot  was  chosen  because 
of  the  extremely  fine  field  date  available  on  it.  The  calculation,  completed  and 
reported  (Ref  1  and  2),  agreed  with  the  field  data  to  within  experimental  err*. 

I  an 4  so  has  given  considerable  confidence  In  these  codes  to  % .jy] ate  atmospheric 

detonations. 

!  te  support  of  DISTANT  PLAIN  pre-shot  planning,  SAP  and  SHELL  were  modified 

to  calculate  the  phenomenology  of  atmospheric  HE  detonations.  A  “burn"  routine 
was  included  in  SAP  to  provide  the  initial  condition#  fo*-  the  calculation  of 
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hydrodynamic  shock  grov  h.  Besides  providing  the  initial  conditions  and  pressure- 
distance  curves,  output  froa  SAP  was  also  used  as  input  to  SHELL,  which  calculated 
tbs  two-dime nalonal  pheaone na  such  as  shock  reflection,  each  stem  formation,  and 
triple-point  path. 

There  were  two  calculations  each  with  a  different  energy  source.  The  first 
was  a  32 -who rt -tom  spherical  charge  of  TXT  detonated  at  a  height  of  26  :-  ters 
Which  was  670  meters  above  sea  level.  The  second  was  a  2 3. 2 -short -ton  spherical 
charge  of  methane-oxygen  mixed  in  a  mole  ratio  of  1  to  1.5  detonated  at  the  aerne 
height  as  the  TXT. 

loth  calculations  were  run  on  SAP  and  SHELL,  using  SAP  output  as  Initial 
conditions  for  the  SHELL  calculation.  The  main  features  of  the  output  are  com¬ 
pared  In  ae  great  detail  aa  possible,  and  the  hydrodynamic  quantities  are  plotted 
for  e  large  number  of  tlatee  In  Volume  II. 
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SECTION  II 
THE  SAP  CODE 

SAP  is  »  one-dimensional  Lagrangian  hydrodynamic  coat.  It  get*  its  name 
(Spherical  Air  PUTT)  from  another  AFWL  one-diocnsio’*  Lagrangian  hydrodynamic 
code,  PUFF  (Ref  3  and  4)  from  which  it  evolved.  FIVE  in  turn  evolved  from  the 
LRL,  Livermore,  SHARP  program.  (Ref  5.) 

SAP  ia  not  a  general  one-dimen* ional  code  bit  ha*  been  specialired  to  prob¬ 
lem*  of  atmoapheric  blast  in  order  to  increase  it*  speed  and  accuracy.  Written 
in  the  spherical  coordinate  system,  it  can  calculate  hydrodynamic  shock*  at  any 
angle  to  the  horizontal.  It  has  a  capability  of  over  100C  rones  art  is  normally 
run  without  sn  automatic  rezone.  The  mesh  is  rezoned  periodically,  ccxsfeining 
zones  "by  hand"  when  necessary. 

Other  features  of  SAP  are: 

a.  Padiation  loss  is  possible,  i.e.,  a  certain  amount  of  energy  may  be 
removed  from  every  zone,  each  cycle,  b}  a  subroutine  that  determines  the  emission 
rate  as  a  function  of  internal  energy  and  density.  (For  the  calculations  reported 
here,  this  option  was  not  uaed.) 

b.  A  "continuous  bun"  i «.  tine  based  on  Chaws*’'- Jouguet  theory  was 

developed  for  SAP.  This  routine  provides  the  conditions  existing  in  the  gaseous 
explosive  products  at  the  ii.rtart  the  detonation  we;*e  reaches  the  surfs  e  of  the 
explosive  charge.  conditions  provide  more  realistic  starting  conditions  for 

the  blast  calculation  than  an  isothermal  sp’..^  .  which  is  more  appropriate  as  a 

starting  condition  for  a  simuiattd  nuclear  burst  calculation. 

c.  Shocks  *-e  calculated  in  the  Von  Keur.nn-Richtmeyer  unm.r.  SAP 
contains  provisions  for  the  use  of  s  linear  and  or  quadratic  a  .ificial  viscoaitv. 

d.  SAP  use*  a  real  atmosphere  that  is  stable  under  *n  8*’  gravity  field, 
and  an  t^uatic®  of  state  for  air  thet  is  an  empirical  fit  to  Hilsenrath's  data 
(Ret  €  and  7). 
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The  equation*  solved  in  SAF  ere  the  partial  differential  equations  fcr 
nonviacous ,  nonconducting  rc&ipressib}e  fluid  flow  written  in  Lagrangian  form. 
These  equations,  including  the  equation  of  state  are; 


Maas 


(*)..*  •  (4  -  * 


(1) 


Momentum 


/l»\  + 1  (h\ 

\9tk  pHt 


(2) 


Snsisi 


I  *A  .  faA 

14 + pl4 


(3) 


Equation  of  State 


p  -  p(p,I) 


(4) 


where  p  -  density  in  gms/cm3 
u  "  velocity  in  cm/sec 
p  •  pressure  in  dynes/cm2 
I  «  specific  internal  energy  in  ergs/gm 
V  ■  specific  volume  in  cm3/gm 
•  1/p 

s  -  energy  source  in  ergs/gm/ sec 

x  ■  Eulerian  coordinate  i..  cm  (position  of  a  point  of  fluid  relative  to 
some  external  reference  frame) 

x  •  lagrangian  coordinate  (permanent  lable  attached  to  each  moving  point 
0  of  fluid) 

g  »  acceleration  due  to  some  external  body  force  in  cm/sec2  (e.g.  gravity) 
t  -  time  in  seconds 


and  where  the  subscripts  above  denote  what  is  being  held  constant  <n  each 
derivative. 


The  finite  difference  approximations  to  the  above  equations,  us  used  in  SAP, 
are  obtained  In  the  usual  manner.  The  fluli  is  divided  into  a  mesh  of  fluid 
e*ew«nts.  Pressure*,  densities,  and  specific  internal  energies  arc  defined  at 
zone  centers.  Velocities  and  positions  are  defined  at  zone  boundaries.  See 

figure  1. 
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The  finite  difference’ equations  are  explicit  time  and  space  centered  equations. 
The  time  centering  is  achieved  in  an  explicit  scheme  by  defining  the  velocity 
and  coordinates  of  the  mesh  points  a  half  cycle  apart  In  time.  The  finite  dif¬ 
ference  equations  are  written: 


Momentum 


(n+1 /2) 


„  to'1/!)  * 

1 

fjn)  (n)  \  / 

lP1-l /2  “  P1+l/2/  4  \ 

(n-1/2)  (n-l/2)\ 

q1-l/2  "  qj+l/2  )  _ 

“3  + 

1/2. 

ft(n)  /  <n)  (n)\ 

Pj-l/2  \  j  J-J 

.  (n)  Mn)  ( n)V 

+  pJ+l/2  V  J+l  "  *j  /. 

Transformation 


At(R+1/2)  +  At(n'1/2) 


x4(n+1)  -  x4(n)  +  U;n+l/2)  At(n+l/2) 

i  j  j 


Continuity 


(n+1) 

Pj-l/2 


(O)  [/  (o)\3  .  /  (o)\ 

.  _ r±±L 


(*ry  ■  (-ct 


+  g 


(5) 


(6) 


(7) 


Znerxz 


r  (n+1  ) 
J-l/2 


n  +  n(n)  +  2n(afl/Z) 

T(n)  _  P2  *  P1-l/2  +  2ql-l/2 

j-l/2  pt  -  p2 

2  ‘  (n) 

P  j  —  1  /2 


(n+1/2)  (n+l/2) 

4  8j-l /2  At 


,(n+l)  v(n)N 

1  "  J 


(8) 
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lauatlon  of  Stats 


“j-l/2  Tj-1/2  h-W) 


m 


where  the  energy  equation  is  the  simultaneous  solution  of  the  first  lew  of  thermo¬ 
dynamics, 

!<*«)  „  j(n)  _  1/2  |p(a+l)  +  p<n^dV  (10) 

and  the  equation  of  state, 


p<n+i)  .  p(»("+i),  i(c+>>) 


and  where 


P2  .  p  („<"+*>.  I<°>) 

p,  -  p  (»(°+1).  i(n)  -  P(n>dv) 

See  reference  4. 

The  artificial  viscosity  term  is  given  by 


pAu 


(—Cj^cs  +  (c0)2  AuJ  if  ff  I  0 


,  3u  _ 
“d  *<0 

«!!>° 

orta2-° 


i.’hich  in  finite  difference  form  becomes 


(11) 

(12) 

(13) 


(14) 


Y  «  '  *  /  * 

Vl/2 


/(n+1/2) 

,  ,(o+l/2)\ 

„  „„(»>  *  .  2  .10+1/2)  (o+l/2)\l 

V'J 

-v>  ) 

-C1  ”j-l/2  +  C0  (“j  -  V,  )J 

if 


(V£/2  -  V£i/l)  i  0  ““  (“I 


(n+1/2) 


.£■'*>)  <  0 


•nd  q5n+l/2)  -  0 

j-l/2 


if 


‘  v“/,)  *  0  «  (“j 


(n+1/2)  _  o(n+i/2)^  >  0 


(15) 


(16) 
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where  cs  is  the  local  sound  speed,  is  the  coefficient  of  linear  viscosity 
(“  0.5),  cQ  is  the  coefficient  of  quadratic  viscosity  («1.8).  The  subscripts 


(j's)  above  are  the  Lagtengian  coordinates  and  the  superscripts  (n's)  refer  to 
time  such  that  t(n+i>  -  t(n)  + 


The  SAP  calculation  advances  explicitly  in  steps  or  cycles.  That  is,  the 
zone  quantities  are  calculated  for  the  next  time  in  terms  of  those  at  the  present 
time.  The  order  of  calculation  proceeds  in  the  same  sequence  as  the  finite  dif¬ 
ference  equations  are  presented: 


a.  A  time  step  is  calculated  such  that  At^0*1^  is  the  minimum  of  all 


!i — hi 


2c  cs  +  2c  2  Au  +  cs 
1  0  1  1 


and  At(n"l/2)  (Ref  3) 


b.  The  momentum  equation  is  solved. 


c.  Transformation  is  accomplished. 

d .  q  ’ 8  are  calculated . 

e.  The  energy  equation  la  solved. 

f.  p^^'s  are  calculated. 


g.  Edits  are  performed  if  desired. 

SAP  being  a  Lagranglan  code  enjoys  all  the  advantages  and  suffers  all  the 
disadvantages  characteristic  of  Lagranglan  codes.  Fine  resolution  is  possible, 
since  those  regions  requiring  -fine  resolution  retain  their  fine  zoning  since  the 
coordinate  system  moves  with  the  fluid.  For  this  reason  fine  resolution  is  found 
in  shock  regions.  However  this  is  at  the  expense  of  resolution  in  rarefaction 
regions.  The  motion  of  material  interfaces  is  also  closely  followed,  since 
coordinate  lines  can  be  placed  along  them  and  will  forever  follow  them. 


Burn 

The  phenomenon  of  detonation  from  the  time  of  its  discovery,  has  been 
theoretically  treated  on  the  basis  of  shock-wave  theory.  Earliest  observations 
of  this  phenomenon  found  that  the  fundamental  property  of  detonations  is  that 
the  detonation  wave,  for  a  given  initial  state  of  explosive  substance  or  mixture, 
is  propagated  through  that  substance  or  mixture  with  a  constant  velocity.  This 
property  has  been  explained  by  Chapman  and  Jocguet  whose  work  has  resulted  in  a 
hydrodynamic  theory  of  detonation. 
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Starting  with  the  experimental  fact  that  datonation  la  a  atable  and  sta¬ 
tionary  process  (balng  suatalnad  by  tha  energy  of  chemical  reaction) ,  they  hava 
ahown  that  tha  datonation  velocity  la  tha  nlnlaun  of  all  possible  velocities  for 
stationary  conditions  of  propagation  of  a  chaaical  reaction  by  a  shock  wave. 

That  state,  corresponding  to  tha  alniaua  value  of  velocity,  has  several  Interest¬ 
ing  properties:  tha  entropy  hoe  a  nlnlaun  value  on  tha  Hugonlot  curve  and  a 
maximum  on  tha  straight  line  In  tha  p-v  plana,  which  joins  the  corresponding 
point  to  the  initial  point,  pQ,  vQ .  The  detonation  velocity  in  this  state  1* 
equal  to  tha  sua  of  the  velocity  of  the  burned  gas  and  the  velocity  of  sound  In 
_£ha-fettrnad-fee . 

Therefore,  baaed  oh  Che  above  properties  of  a  Chapman-Jouguet  detonation, 
a  "burn"  routine  was  wrlttan^for  SAP.  Two  methods  of  "burn"  were  used,  both  of 
which  gave  excellent  results.  In  the  first,  the  theoretical  detonation  velocity 
and  the  energy  released  per  gran  of  explosive  substance  are  entered  as  Input 
nuabera.  Then 

„<«•■)  .  ,<n>  +  Dit(»+l/2)  (17) 


where  R  »  radius  o S  the  detonation  wave  In  an 
D  -  detonation  velocity  In  ca/sec 
At  »  time  step  in  sec 

AI  -  energv  released  during  the  time  step  In  ergs 
p  ■  density  of  the  cell  which  Is  currently  being  "burned”  in  gms/cm3 
k  -  Index  of  cell  currently  being  "burned" 

Q  "  energy  liberated  during  detonation  In  ergs/gm 

In  the  second  method  only  the  energy  released  per  gram  of  explosive  substanc. 
detonated  Is  entered  *s  an  input  nunber.  When  the  transformation  calculation  is 
performed  the  radius  of  the  detonation  wave  is  similarly  jpdated: 


R(n>  + 


(n+l/2)  (n+1/2) 

>+■  /*>  _ iVl _ 

\-l  (n+l) (nf 

\  '  Vl 


.(n+l) 


At(n+l/2)  (19) 


Then  the  new  radius  of  the  detonation  wave  and  the  energy  released  Is 
calculated: 


• 


AFWl-TR-66-141,  Vol.  I 


-  i  +  c.At(B+1'S> 


41  -  ^R("+J,)J  -  it3j 


(n+l) 


(20) 

(21) 


where  cs  la  the  sound  speed  In  the  burned  gaa. 

Therefore,  the  sequence  of  calculation  with  "burn"  becomes: 


Calculate  t 


(n+l/2) 


(n+1/2) 

J 


b.  Solve  the  momentum  equation  for  all  u 

c.  Update  the  Eulerlan  coordinates:  ■  xjn^  +  yfo*1/2)  t(n+l/2) 

d.  Calculate  R  aa  shown  above  If  method  2  "burn"  la  used. 

e.  Calculate  the  new  radius  of  detonation  wave. 


Calculate  energy  released  due  to  the  detonation  of  material  In  this 
time  cycle  and  add  the  energy  to  the  cell. 


g.  Sclve  the  energy  equation. 

h.  Calculate  the  new  pressures  from  the  equation  of  state. 

1.  Go  to  a  and  begin  next  cycle. 

This  sequence  Is  followed  until  all  the  explosive  matter  h*»  been  detonated. 
Then  the  calculation  becomes  a  "pure"  hydrodynamic  calculation. 


4 
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SECTION  ZIX 
THE  SHELL  CODE 

SHELL  It  t  versatile  two-dimensional  hydrodynamic  code.  It  can  b«  used  to 
solve  pxoblem*  In  any  medltat  whose  notions  can  ba  analysed  on  the  baele  of 
invitcid  hydrodynamic  behavior  of  a  cyllndrlcally  symmetric,  conpreaelble  fluid. 
Some  typical  problem  ere:  hypervelocity  impact,  cratering,  vehicle  re-entry, 
and  detonations,  underground,  underwater,  and  in  the  atmosphere.  A  model  and  an 
equation  of  state  for  the  aediua  of  the  problem  are  required  by  the  code. 

The  Air  Force  Weapons  Laboratory  has  two  versions  of  the  SHELL  code.  Both 
were  developed  at  General  Atomic  by  b.  E.  Freeman  and  W.  E.  Johnson  primarily 
for  hype  oclty  Impact  calculations.  The  codes  were  modified  at  AFWL  to 
calculate  rircball  rise  and  expansion. 

The  firat  and  older  of  the  two,  SHELL-PIC,  uses  the  Partlcle-in-Cell  Method, 
developed  by  F.  H.  Harlow  (Ref  8)  at  Los  Alamos  In  1955,  to  solve  the  hydro- 
dynamic  equations.  ’  The  second,  SHELL-OIL,  Is  a  continuous  version  of  SHELL-PIC. 
Since  early  1963,  development  of  SHELL-OIL  has  been  continued  by  J.  M.  Walsh  and 
W.  E.  Johnson. 

SHELL-PIC  and  SHELL-OIL  have  been  in  use  at  AFWL  since  1961  and  1964, 
respectively. 

In  the  following  pages  only  SHELL-OIL  will  be  described,  since  it  was  used  in 
the  calculations  reported  here. 

The  equations  solved  in  SHELL  are  the  partial  differential  equations  for 
nonvlscous,  nonconducting  compressible  fluid  flow.  These  equations,  including 
the  equation  of  state  ar~ 

Mass 


+  0  ♦  pV-u  *  0 

(22) 

Momentum 

p(j~  +  1j  ♦  Vp  +  pV*  -  o 

(23) 

Eneray 

°(n  +**v)  E  4  m  0 

(24) 

Equation  of  State 

p  -  p(?,i) 

(25) 
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where  p  =  density  in  gm*/cm3 
^  =  velocity  in  ca/sec 
p  =  pressure  in  dynes /cm2 
E  =  specific  total  energy  in  ergs/ga 
-1+1/2  u2 

I  =  specific  Internal  energy  in  ergs/ga 
t  =  tiae  in  sec 

t  s  potential  of  external  force  field  In  ergs/ga 
a.  SHELL-OIL 


SHELL-OIL  is  a  one  aaterial ,  pure  Euler lan  code.  It  solves  the  hydro- 
dynamic  equations  by  dividing  the  region  occupied  by  the  fluid  Into  a  aesh  of 
fixed  cells.  The  fluid  is  then  described  at  any  instant  of  tlae  by  specifying  the 
velocity,  density,  and  internal  et  xgy  for  each  cell.  These  values  are  considered 
to  be  known  at  the  center  of  each  cell  end  constant  over  it.  This  forced  hoao- 
genelty  of  the  cell  introduces  a  false  diffusion,  one  disadvantage  of  Euler lan  codes. 

Other  disadvantages  of  SHELL-OIL  are  the  general  inability  of  Euler las 
codes  to  resolve  fine  detail  Moving  with  the  fluid,  since  the  aesh  is  fixed  in 
space  and  the  Inability  to  follow  aaterial  interfaces.  However,  Euler lan  codes 
enjoy  the  major  advantage  of  being  able  to  solve  the  bydrodynaalc  equations  even 
in  the  presence  of  large  fluid  distortions. 

To  begin  a  SHELL-OIL  calculation,  the  problea  must  first  be  generated  by 
CLAM,  the  generator  code  for  SHELL  .  CLAM  sets  up  the  aesh,  giving  each  cell  disten¬ 
sions  and  the  following  quantities:  a  radial  and  an  axial  velocity  component,  a 
mass,  and  a  specific  internal  energy. 

The  cell  quantities  are  derived  froa  data  entered  in  groups  of  packages. 
These  data  include  the  type  of  eaterial,  dimensions,  velocity  components,  density, 
and  specific  internal  energy  of  the  package.  To  describe  as  aany  geometries  and 
energy,  density  and  velocity  distributions  as  possible,  CLAM  places  N2  particles 
into  each  cell,  where  1  <_  N  <  20  end  is  also  included  in  the  package  data.  Each 
cell  ia  divided  into  N2  equal  parts,  end  e  particle  ie  placed  at  the  center  of  each 
area.  Taking  this  package  data,  a  fit  assigns  each  particle  two  velocity  com¬ 
ponents,  a  density  and  a  specific  Internal  energy.  Then  the  naas  of  each  particle 

ie  computed:  the  density  tiaes  the  volume  of  that  part  of  the  cell  containing 
the  particle.  The  mas*  of  the  cell  ie  the  sun  of  the  aaasce  of  ell  the  particle# 

in  the  cell.  Soth  nomen tun  components  are  calculated  as  the  sum  of  the  lnifvldual 

momentum  components  of  each  particle  in  the  cell.  The  internal  energy  of  each 
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Mil  is  Um  sun  of  tho  Internal  energies  of  all  the  particles  in  the  cell. 
Finally,  tbeee  call  quantities  are  converted  to  two  velocity  components  and 
•pacific  Internal  energy. 

the  goo— trlae  available  la  CLAM  arat  rectangle,  triangle,  ellipse,  per¬ 
turbed  ellipse  and  circle.  The—  geo— tries  are  In  the  r-s  plane  and  take  their 
3~di— ional  counterparts  upon  rotation  about  the  z-axls. 

After  all  paebagee  have  bean  procaaaad,  CLAM  puts  all  call  dl— Iona  and 
quantities  on  tap#.  (If  CLAM  were  generating  the  problen  for  a  FHELL-PIC  calcu¬ 
lation,  the  coord  loot ee  and  masses  of  each  particle  would  also  be  written  on  thla 
tape.)  This  output  provides  the  starting  conditio—  for  tho  SHELL  calculation. 

Since  SHELL-OIL  la  e  on*  notarial  code,  an  option  to  inaert  — aelaaa  trace 
pert  Idea  into  the  problen  has  bean  included  in  CLAM  at  AFVL.  The  purpose  of 
the—  trace  pert  idee  ie  to  follow  the  detailed  notion  of  e  labeled  —  terlal,  to 
el— late  the  notion  of  e  Lagranglen  interface,  or  both.  Therefore,  the  CLAM 
output  also  include*  the  coordinates  of  these  trace  pertidee  on  tape  in  addition 
to  the  call  quant  it  lee  end  di— ions  for  e  SHELL-OIL  calculation.  These  trace 
particles  do  — t  affect  the  hydrodynamics  in  any  umur. 

The  SHELL  calculation  advances  explicitly  in  steps  or  cycles.  That  is,  the 
cell  quantities  are  calculated  for  tine  t  +  At  in  ter—  of  chose  at  tine  t.  The 
solution  of  the  hydrcdynaaic  equations  is  divided  Into  two  phases.  In  the  first, 
the  Euler lan  f.sld  functions  ere  updated  considering  pressure  effects  slone, 
while  neglecting  —terlal  transport  due  to  fluid  notion.  In  the  second  phase, 
continuous  mass  flow  across  cell  boundsrlbs,  new  c*.'l  densities,  velocities,  and 
specific  internal  energies  are  calculated.  Hess,  no— n tun,  and  total  energy  ere 
conserved .  If  —as  flows  out  of  the  —eh,  an  nuto—  tic  rezooe  increases  the  size 
of  the  —eh. 

There  ere  three  tine  control  conditio—  thst  govern  the  cl—  step,  A t.  They 

are: 

e.  The  coursnt  condition 

b.  The  asxinui  of  IjJfj  and  !^l  « 

c.  A  negative  tine  step  opt  ion 

The  first  two  conditions  prohibit  the  tranan lesion  of  a  eignel  or  nsss  ecroe- 
no re  then  cne  cell  In  one  tine  interval.  This  require— nt  ie  imposed  ly  con¬ 
sideration  *  of  stability  of  the  finite  difference  equations  (kef  8).  The  third 
condition  will  be  discussed  below. 
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After  a  time  step  ha*  been  determined,  SHELL  perform*  the  first  phase  of  the 
hydrodynamic  calculation  where  the  Eulerian  description  of  the  fluid  Is  advanced. 
This  is  done  according  to  a  finite  difference  approximation  to  the  hydrodynamic 
equations,  Eq  22,  23,  and  24 

Conservation  of  mass  (Eq  22)  Is  automat l'** l ly  satisfied  by  the  fluid  model. 
Hass,  which  leaves  one  cell,  enters  another  and  the  change  in  mass  is  accordingly 
subtracted  from  the  donor  cell  and  added  to  the  receiver  cell.  Mass  is  neither 
created  nor  destroyed. 

Conservation  of  momentum  and  energy,  (Eq  23  and  24),  ia  treated  as  follows: 
the  first  phsse  considers  the  fluid  st  rest  snd  determines  only  the  contribution 
of  the  pressure  terms  to  the  time  derivative.  Therefore,  the  transport  term 
(second  term  on  the  left-hand  aide)  of  each  equation  is  drooped.  Since  the 
hydrodynamic  quantities  are  defined  at  the  cell  centers,  the  equations  in 
finite  difference  form  become: 


Momentum 


where 
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(26) 


(27) 


(28) 


and  k  ia  the  index  of  tbe  cell  center  and  i  la  the  index  of  the  right-band 
boundary  of  cell  k.  See  figure  2. 

Two  passes,  in  succession,  are  made  through  tbe  first  phase  of  calculation  of 

each  cycle.  In  tbe  firet  pass,  tentative  new  vslues  of  velocity  sre. calculated 
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for  a  fall  t  1m  step.  The  Internal  energy  la  updated  by  half  a  tine  atop,  using 
t he  old  velocities.  Is  the  second  paas ,  tha  lntaraal  energy  la  updated  another 
half-tine  step,  using  the  tentative  new  velocities  calculated  la  tha  first  pass. 
The  a&aerey  eduction  la  treated  In  this  fashion,  as  proposed  by  Harlow  (Xaf  8) , 
feats  considerations  of  stability  as  wall  aa  greatest  accuracy  in  the  behavior  of 
fluid  an. ropy.  Tha  anargy  of  tha  systan  la  conserved  Identically,  even  in  the 
finite  difference  approach. 


In  both  pa  seen,  tbs  lntaraal  energy  la  checked  to  aaa  If  the  updated  value 
baa  became  negative.  If  it  hae,  an  option  exists  to  integrate  back  to  the 
original  configuration  of  tha  fluid  at  tha  beginning  of  tha  cycle  by  using  the 
negative  At.  Than  the  Integration  goes  forward  again  with  a  mailer  At  which 
has  bean  calculated  add  which  will  knap  the  internal  energy  positive, 
d  mot 
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*tc'  Figure  2.  The  SHELL  Mesh. 

Ia  tha  secosd  phase  of  the  hydrodynamic  calculation  ccrattnaou*  mass  trans¬ 
port  is  calculated  according  to  the  equation  for  the  conservevioo  on  mass,  which 
in  f'  *  difference  form  ia: 


A£  Pi-iViVi 

tl  ri-i/lir 


rlclui 
*1-1 /2*1 
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This  can  be  rewritten  to  the  following  form,  which  is  used  by  the  code. 


Mkn+0  *  "I!10  +  VLt  [('i  "  '  Vt)+  2AZ  fi-lVl 

Where  v^p^  i  mass  flu*  across  bottom  face  of  cell  K 
-  mass  flux  across  top  face  of  cell  K 
ulpl  J  mass  flux  across  left  face  of  ceil  K 
u^Py  =  mass  flux  across  right  face  of  cell  K 


riV* 


)] 


(30, 


and  where  the  densities  are  t!  -se  of  toe  donor  cells  and  the  velocities  are 

I 

1/2  (v  +  \) 
v'~ 

1/:K +  \) 

'mTTfar%* 

where  M  signifies  the  index  of  the  donor  ceil.  The  velocity  weighing  scheme  used 
above  is  one  that  en  ures  stability  in  regions  behind  a  shock  front  (Ref  9). 


When  mass  transport  is  being  calculated  for  ?  given  cell,  the  transfer  at 
the  left  and  bottom  faces  is  already  available  from  the  previous  column  sweep 
and  cell  lolow,  respectively.  Therefore,  for  a  particular  cell,  the  mass  trans~ 
port  is  calculated  for  the  top  ana  right  feces  only.  After  the  masses  trar.  - 
ferred  at  all  faces  of  a  cell  have  been  determined ,  the  a  went a  and  the  total 
srttifi''  energy  associated  with  these  masses  are  corpured.  Then,  by  conserving 
both  axial  and  radial  momentum,  the  new  veloclti  ire  calculated.  The  *c?tal 
energy  carried  by  the  transported  masses  is  subtracted  from  U:^  donor  cell  and 
added  :o  the  receiver  ceil.  Then,  by  conserving  total  energy,  the  new-  internal 
energy  of  each  ccJi  is  the  dfilerenee  between  T.  he  n?v  total  ei.irgv  anti  the  new 
kinetic  energy.  The  new  specific  Internal  energy  is  the  new  internal  rwr!'^ 
•Jiviaed  by  the  nev  mass. 
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t(n+l 

*k 


’  •  |Mk">[1k  +  1/2(JU  +  Jk)j  +  4Mb[*M  +  l/2('iM  +  '4 

+  “l^M  +  lri^y  *  »«)] '  +  H  +  s)] 


(3«) 


I 


Whet  £  AMg,  AHj,»  end  are  the  masses,  in  gms  transferred  across  the  bottom, 

left,  top,  and  right  faces  of  cell  K,  respectively,  and  where  the  quantities 

enclosed  in  brackets  and  subscripted  by  M  are  he  specific  total  energies  of  the 

donor  cells.  The  above  are  the  flral  values  of  velocity,  mass,  and  specific 

internal  energy  for  the  cycle.  , 

if  any  mas:  leaves  the  top,  right,  or  bottom*  boundaries,  thac  mass  and  the 
associated  energy  are  subtracted  fro-i  the  total  mass  and  energy  of  the  system 
respectively.  If  the  mass  leaving  the  r^osh  is  a  certain  percentage  of  the  mass 
of  the  cell  from  which  it  left,  a  flag  is  set  to  rescue  the  system.  Additional 
flags  are  also  set  to  indicate  whether  the  top,  bottom,  and  rignt  boundaries, 
or  any  combination  of  these  three,  of  Che  mesh  have  been  crossed. 

The  automatic  rezonr  for  SHELI -GIL  has  been  developed  at  AFWL.  In  tezone, 
a  predetermined  number  of  columns  and  or  rows,  depending  whether  the  top,  bottom, 

I 

and  or  right  boundaries  have  been  crossed,  are  deleted.  The  remaining  ceils  are 

expanded  in  else  to  fill  the  original  geometry  of  the  mesh.  The  mass  or  each 

expanded  cell  is  the  sum  of  tne  masses  of  those  portions  of  the  original  cells 

which  raeke  up  the  new  cell.  By  conserving  momentum  and  total  energy,  new 

velocities  and  internal  energies  of  the  new  ceils  are  calculated.  Than  the  same 

number  of  columns,  rows,  or  both,  that  ha^e  previously  been  deleted  are  now  added 

to  the  mesh  thus  regaining  the  original  number  of  ceils.  Ambient  conditions  are 

give,.  to  these  added  cells  The  calculation  ia  then  resumed.  I 


The  calculation  of  mass  transport  ir  SHELL-OIL  has  anr.her  feature  that 
removes  preferential  transfer  caused  by  initial  choice  of  indexing.  If  the  mass 


There  Is  _.lso  an  option  for  a  reflective  bottom  boundary  as  was  use^  in  the 
aK.jlati.on  reported  here. 
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out  the  top  and  right  would  remove  more  than  the  mass  in  the  cell,  the  code 
recalculates  new  mess  transfer  by  a  weighing  procedure.  The  AMg  would  be  Its 
fraction  of  the  total  mass  out  times  the  mas «  of  the  cell  and  AHj,  would  be  its 
fraction  of  the  total  mass  out  times  the  mass  of  the  cell.  That  Is, 


smt  '  l^hir  (4n)  +  s\  *  smb) 

aMR  •  is^r3£  (4n)  +  +  Js) 


(35) 

(36) 


In  this  calculation  AM^  and  AMg  are  considered  zero  if  they  are  flowing  into  cell 
K  so  thar  the  maximum  mass  that  could  flow  out  the  top  and  right  is  .  However 
AM^  and  AMg  are  flowing  out  of  cell  K,  the  cell  la  committed  to  deliver  this 
mass  and  so  the  maximum  mass  that  can  flow  out  the  top  and  right  is  M^n^  +  AM^ 

+  AMg  (where  and  AMg  are  negative  signifying  flow  ir.  negative  direction). 


After  the  final  velocities  for  the  cycle  have  been  determined  for  all  the 
cells  in  the  mass  transport  phase  of  the  calculation,  trace  particle  movement 
is  calculated.  The  trace  particles  play  no  part  in  the  hydrodynamics  of  the 
system  but  are  merely  moved  to  new  positions  according  to 


Xg  +  U 'At 

(37) 

y0  +  v'At 

(38) 

The  velocities  used  to  move  each  particle  are  the  local  velocities  at  the  location 
of  the  trace  particle.  They  have  values  which  are  Interpolated  between  the  velo¬ 
cities  of  the  cells  whose  centers  surround  the  partible.  The  velocities  are 
obtained  by  assuming  a  rectangle  with  ths  same  dimensions  as  a  cell  to  be  located 
around  the  particle.  This  imaginary  cell  overlaps  from  one  to  four  mesh  cells, 
depending  on  the  position  of  the  parti “le  in  the  mesh  cell.  Then  the  interpolated 
velocity  is  calculated  as  the  weighted  average  of  the  cell  velocities,  the  weigh¬ 
ing  being  proportional  to  the  overlap  areas.  Hence  the  motion  of  the  trace 
particles  follows  the  motion  .f  material,  or  a  Lagrangian  Interface,  depending 
on  how  the  particles  were  originally  placed  in  the  problem,  and  gives  a  picture 
of  the  flov  pattern.  After  all  trace  particjes  have  beer,  moved,  the  calculations 
for  the  cycle  have  been  completed,  and  the  next  cycle  begins  with  the  determina¬ 
tion  of  a  new  time  step. 
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b.  Stability 

Tht  Eulerian  aquations  uaad  ia  tha  firat  phase  of  the  calculation  are 
unstable,  since  they  do  not  contain  tha  dissipative  mechanism  necessary  for  a 
finite  difference  technique  to  calculate  shocks.  Since  shocks  appear  as  mathe¬ 
matical  surfaces  on  which  aueb  fluid  properties  as  density,  pressure,  Internal 
energy,  and  entropy  have  discontinuities,  suitable  boundary  conditions  (those  pro¬ 
vided  by  the  Sanklne-Sugonlot  equations)  are  needed  to  connect  values  of  th» 
above  quantities  on  both  aides  of  the  shock.  However,  Von  Neumann  and  Richtmeyer, 
(Ref  10),  have  shown  that  the  hydrodynamic  equations  can  be  straightforwardly 
solved  by  numerical  methods  if  an  artificial  dissipative  term  is  introduced. 

This  tens  smears  the  ehock  so  that  the  surface  of  discontinuity  is  replaced  by  a 
thin  layer  in  which  the  above  quantities  vary  rapidly  but  continuously.  Hence, 
the  numerical  calculation  proceeds  as  if  there  were  no  shock  at  all  and  at  the 
same  time  satisfying  the  Kankine-Hugcnlot  conditions. 

The  instability  of  the  calculations  in  the  first  phase  introduces  the  main 
source  of  entropy  loss  in  the  entire  calculation.  However,  stable  cflculations 
can  be  made  by  S3£LL,  since  errors  introduced  in  the  second  phase  have  the  effect 
of  smoothing  out  the  discontinuities. 

Harlow  (Ref  8)  has  shown  that  the  treatment  of  mass  movement  in  the  second 
phase  of  the  calculation  produces  dissipative  effects  that  give  stability  to  the 
calculation.  He  demonstrated  this  by  expanding  the  difference  equations  in  a 
Taylor  series  about  some  central  space  and  time.  The  result  was  the  original 
differential  equations  plus  some  additional  terms  thst,  in  lowest  order,  have 
the  appearance  of  true  viscous  and  heat  corsductijn  terms.  Hence  th-sy  are  called 
the  effective  viscosity  cud  effective  heat  conduction. 


Conservation  of  Mass 


+  llrpu)  3(py) 
3t  r3r 


Conservation  of  Momentum —  direction 
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Conservation  of  Momentum — z  direction 


/3v  3v  3v\  3p  3  /  .  3v\  /  3v  \ 

p\3?  +  “  Ti*  v  Til*  fi  m  757  lrX„  Ti)*Ti  1‘,  Til 


to* 


AFWL-TR-66-141,  Vol.  I 


t 


Conservation  of  Energy 

+  \  [(It)2  -  m3]  (42) 

where 

*u  •  l/2p |u|&r 
*v  -  l/2p |v|  Az 

However,  Bjork  (Ref  11)  has  shown  tLat  any  attempt  to  relate  the  terms 
appearing  on  the  right-hand  side  of  the  equations  above  to  distinct  physical 
effects  leads  to  many  contradictions.  Some  are:  the  "effective  viscosity"  terms 
contain  symmetric  and  antisymmetric  elements ;  the  c  igin  of  the  symmetric  elements 
is  not  viscous  in  nature;  the  presence  of  the  antisymmetric  elements  just  adds 
confusion  to  any  attempted  representation;  the  terms  do  not  leave  the  flow 
equations  in  an  invariant  form;  a  heat  conduction  law  is  implied  where  heat  flow 
is  proportional  to  the  gradient  of  specific  internal  energy  rather  than  of 
temperature.  Moreover,  the  values  of  these  terms  far  exceed  the  magnitude  of 
the  analagous  physical  terms.  Therefore,  these  terms  cannot  be  thought  of  as 
anything  except  errors,  but  they  do  make  SHELL  a  practical  code. 

The  form  of  the  effective  viscosity  suggests  another  limitation  of  the  SHELL 
code.  Since  the  viscous  effects  are  proportional  to  the  local  mean  fluid  speed, 
in  regions  of  low  fluid  speed  the  viscosity  will  become  ineffective  and  the 
instabilities  of  the  difference  equation  will  cause  an  exponential  growth  of  any 
perturbations  to  the  solution.  But  as  the  velocities  increase  and  become  com¬ 
parable  to  the  local  sound  speed,  the  viscosity  again  takes  effect  to  limit 
further  growth  of  instability.  Hence,  the  instability  is  bounded  and  calculations 
can  be  made.  This  effective  viscosity  is  the  main  source  of  entropy  production 
in  the  calculation. 
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SECTION  IV 
THE  CALCULATION 

The  calculations  reported  here  were  performed  on  the  CDC-6600  computer.  Cade 
development,  calculation,  and  output  analysis  (plotting)  required  approximately 
100  hours  of  computer  time. 

1.  SAP-TNT 

The  SAP-TNT  calculation  was  taken  from  0  to  0.534  second.  Initially,  800 
cones  were  used.  The  first  cone  had  a  Ar  of  0.12  cm,  and  the  other  799  had  Ar 
of  0.5  cm.  The  first  327  cones  were  TNT  (32  short  tons*) .  The  initial  condi¬ 
tions  were:  density — 1.56  gm/cm3,  Internal  energy — 0,  velocity — 0.  Beyond  these 
327  cones,  the  remaining  cones  were  the  ambient  atmosphere  (670  meters  altitude). 

"Burn"  method  1,  as  described  previously,  was  used  to  detonate  the  TNT.  The 
detonation  velocity  (corresponding  to  a  loading  density  of  1.56  gms/cm3)  used 
was  6.81  *  10s  cm/sec.  The  energy  released  per  gram  of  TNT  detonated  was 
4.264  *  1010  ergs/gm  (Ref  12). 

The  LSZK  equation  of  state  for  the  detonation  products  of  TNT  (Ref  12)  was 
used  for  the  TNT  cones  that  were  "burned."  For  the  unburned  TNT,  a  compressions! 
equation  of  state,  which  ls«an  analytical  fit  to  the  Hugonlot  for  solid  explo¬ 
sives,  was  used. 

2.  TNT- SHELL 

SAP  results  at  6.62  milliseconds  were  scaled  to  20  short  tons  of  TNT  and 
used  as  input  for  the  SHELL  caxculation.  As  a  result  of  scaling,  the  SHELL 
starting  time  became  5.66  milliseconds  The  SAP  values  for  density,  specific 
Internal  energy  and  velocity  were  fitted  into  the  two-dimensional  SHELL  mesh 
containing  a  real  (1962  standard)  atmosphere. 

A  mesh  156  cells  in  the  c-direction  and  142  cells  in  the  r-direction  for  a 
total  of  22,152  cells  was  established.  All  cells  except  the  bottom  row  of  142 
cells  comprised  the  active  mesh.  Each  cell  in  the  bottom  row  was  644  meters 
high  and  30  cm  vide.  Each  cell  in  the  active  mesh  was  a  30-cm  square. 

*Thls  particular  charge  was  chosen  to  save  calculation  time,  since  it  was 
directly  applicable  to  a  proposed  surface  detonation  of  a  hemispherical  charge 
of  20  tons  of  TNT  (assuming  20  percent  energy  loss  to  the  ground)  and  scalable 
to  Distant  Plain  event  one. 
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Therefore,  the  overall  mesh  extended  from  se*-l:vel  to  an  altitude  of  690.5 
meters  and  out  to  a  radius  of  42.6  asters;  ar.d  the  active  portion  of  the  aesh 
froa  644  aeters  (ground-level)  to  an  altitude  of  690.5  aeterc  and  out  to  s 
radius  of  42.6  meters.  The  bottom  boundary  (ground-level  at  644  aeters)  was 
reflective;  whereas,  the  top  and  right  boundaries  were  transaittlve  and  thus 
subject  to  rezone  In  their  respective  direction*.  The  left  boundary  (r  ■  0) 
was  the  axis  of  symmetry. 

Data  were  entered  into  CLAM  in  t.fo  packages,  each  of  which  generated  25 
particles  per  cell  for  the  sole  puroonc  of  fitting  the  data  into  the  SHELL  mesh. 
The  first  package  fitted  the  cablent  atmosphere  into  the  active  mesh  with  the 
exception  of  those  cells  conialnec.  in  a  sphere  whose  center  was  on  the  x-axls 
at  an  altitude  of  670  meters  and  whose  ranlus  was  16.55  meters.  Into  these  cells 
the  second  package  fitted  the  S.JP  results  of  velocity,  density,  end  Internal 
energy  at  6.62  aillisectads  (5.66  as  seeled). 

Trace  particles  we -e  introduced  into  the  calculation  in  order  to  follow  the 
TNT-air  boundary.  A  .agranglan  interface  of  921  particles  was  place!  In  a  semi¬ 
circle  to  simulate  ne  TNT-air  boundary.  The  center  of  this  circle  was  on  the 
z-axls  at  an  altitude  of  670  asters,  and  the  radius  was  14.188  meters. 

Also  Introduced  into  the  calculation  were  19  "teat  stations'*  where  the  over¬ 
pressures,  dynamic  pressures,  impulses,  and  velocities  ware  calculated  ur  each 
cycle.  The  positions  of  these  "test  stations"  were  chosen  to  coincide  with 
experimental  stations.  See  table  I  and  figure  3. 

The  equation  of  state  used  in  the  SHELL  calculation  was  that  of  air,  tba  asms 
equation  of  state  for  air  as  used  in  SAP.  This  equation  of  state  was  used  for 
ell  regions  Including  those  regions  that  contained  TNT. 

3.  SAP-Methane 

The  SAP -methane  calculation  was  taken  from  0  to  1  second.  Initially  800 
zones  were  used  in  the  calculation.  The  first  rone  had  a  Ar  of  6.4  cm,  and  the 
other  799  had  Ar  of  5.0  cs.  The  first  335  zones  contained  the  methane  balloon 
(Initial  radius:  16.764  aeters).  The  initial  conditions  were:  density — 

1.07197  *  10~3  gma/ca3,  internal  energy — 2.1886  *  10*  ergs/gm,  velocity— 0. 

Beyond  these  335  zones,  the  remaining  tones  were  the  ambient  atmosphere  (670 
aeters  altitude). 
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"•urn"  method  2,  as  described  previously,  ms  used  to  detonete  the  methane. 
The  energy  releeeed  per  gren  of  aathana-oxygen  mixture  detonated  was  8.01  *  1010 
ergs/gn  (Kef  13).  The  detonation  velocity  calculated  using  "burn"  aethod  2  was 
2.M3  *  10s  ca/sec,  as  compared  to  a  theoretical  2.752  *  10s  ca/sec  (Ref  13). 

The  Ideal  gas  equation  of  state  was  used  for  the  burned  and  unburned  aethane. 
For  burned  aethane,  the  value  of  gaaaa  used  was  1.2136,  while  for  unburned  aethane, 
the  value  of  gaaaa  used  was  1.3986  (Ref  13). 

4.  Met-henm-SHRLL 

SAP -me thane  results  at  8  ailliseconds  were  used  as  Input  to  SHELL.  The  SAP 
values  for  density,  specific  Internal  energy,  and  velocity  as  a  function  of  posi¬ 
tion  were  fitted  into  the  two-dimensional  SHELL  aesh.  The  geoaetry  and  boundary 
conditions  were  the  sane  for  the  aethane  calculation  as  for  the  TNT  calculation 
previously  described. 

Data  were  entered  into  CLAM  in  two  packages,  each  of  which  generated  25 
particles  per  cell  for  the  sole  purpose  of  fitting  the  data  into  the  SHELL  aesh. 
The  first  package  fitted  the  aabient  ataosphere  into  the  active  aesh  with  the 
exception  of  those  cells  contained  in  a  sphere  whose  center  was  on  the  z-axis  at 
an  altitude  of  6/0  aeters  and  whose  radius  was  20.21  aeters.  Into  these  cells, 
the  second  package  fitted  the  SAP  results  of  velocity,  density  and  Internal  energy 
as  a  function  of  position  a*-  8  ailliseconds. 

Again  trace  particles  were  introduced  into  the  calculation  to  represent  the 
aethane- air  boundary.  This  Lagrangian  interface  consisted  of  921  particles  placed 
In  a  aealclrcle  whose  center  was  on  the  z-axis  at  an  altitude  of  670  aeters  and 
whose  radius  was  17.0  aeters. 

"Test  stations"  were  introduced  in  the  aethane  calculation.  The  location  of 
these  19  stations  were  the  saae  ss  for  the  TNT  run.  Again  overpressures,  dynaaic 
pressures,  impulses  and  velocities  were  calculated  each  cycle  *t  these  points. 

See  table  I. 

5.  LSZK  Equation  of  State  (Ref  12  and  14) 

Landau,  Stanyukovich,  Zeldovich,  and  Kompaneets  have  developed  an  equation 

of  state /for  the  explosion  products  of  s  condensed  explosive.  They  treated  the 

* 

highly  compressed  gas  as  s  solid,  likening  the  state  of  the  detonation  products 
immediately  following  completion  of  reaction,  to  the  crystal  lattice  of  the 
•olid  state.  Further  details  concerning  this  equaticn  of  state  are  given  in  the 
above  aentlrned  references. 
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The  LSZK  equation  of  state  can  be  vi It ten 

where  p  ■  pressure  *n  dyne s/c*2 

I  ■  specific  internal  energy  in  erga/ga 

a,v  -  dimensionless  constants  evaluated 
using  experimental  results 

B  »  constant  with  dimensions  (ML-3)  1  -  v  also  to  be  evaluated  from 
experimental  results. 

Lutzky  (Ref  12)  gives  the  values  of  the  above  constants  for  THT  as 
-  2.9412 

B/Q  -  0.53562  1  -  v 

v  •  2.78 

where  0  -  1018  calories /gn  is  the  heat  of  detonation  and  the  ideal  gas  value  for 
the  ratio  of  specific  heats  of  the  products  at  low  densities  was  taken  to  be 
Yj  -  1.34. 

This  equation  of  state,  as  used  in  SAP  for  THT,  becomes 

p  -  0.34  Ip  +  1.877  *  10*V’79 
6.  Equation  of  State  for  Air 

The  SAP  end  SHELL  equation  of  state  for  air  is  a  fit  derived  by  APVL  (lef  7) 
f roe  the  data  of  reference  6.  Since  L.ie  codes  use  energy  end  density  as  the 
independent  variables  and  are  pure  hydrocodea,  the  only  requirement  on  the  equa¬ 
tion  of  state  is  pressure  as  a  function  of  energy  and  density.  It  is  ccwston  to 
consider  an  effective  gaana,  which  ii  1  +  P/p I.  The  equation  cf  state  returns 

-  1)  from  which  is  onstructed  the  pressure.  The  region  of  interest  of  SAP 
and  SHELL  is  chiefly  molecular  air  and  therefore  the  detailed  equation  of  state 
la  q,  te  complicated.  The  present  equation  is  a  fairly  complex  formulation  that 
gives  all  of  the  structure  of  this  molecular  region  although  not  always  with 
extrtrsly  high  accuracy.  The  portion  of  the  equation  of  state  of  interest  in  this 
problem  ia  the  10-2  to  10*  density  region  below  2  '  10l 1  ergs/ga.  In  this  region 
the  fit  la  almost  always  within  5  percent  with  an  average  error  of  about  2  or  3 
percent.  It  is  however  an  interim  equation  an*'  will  he  i ’placed  soon  by  a  more 
accurate  formulation.  The  analytic  expression  of  this  equation  of  state  is  given 
belcw.  I  1*  specific  energy  in  units  of  Jerks /wegag ram  or  10‘°  ergs/g«. 
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7 .  Atmospheric  Hcdel 

The  SH'  LJL-Oli.  er-oepberic  c  Mlai  v as  also  derived  by  AIVL  from  data  contained 
in  reference  15.  Tc  perform  hydrcdj  '*s>ic  calculations  in  thv  atnespbere,  the 
atmosphere  mi  at  be  very  stable,  otberv>*«  tha  firsw  phaae  of  the  calculation  will 
result  in  fictional  vertical  velocities.  Tr.S  *~nditioo  la  given  by  the  hydro¬ 
static  equation 

*  -8**) 


The  atmosphere  is  divided  into  22  Altitude  group*  frees  0  to  7i>o  ka  uhere  the 
molecular  temperature  varies,  approximately  linearly  with  altitude  within  each 
group.  That  i# 


T  £1)  -  T  ♦  L 
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wh^re  th*  subzero  values  refer  to  the  oase  altitude  of  each  group, 
the  following  relationships: 


and 
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the  hydrostatic  equation  becomes 

d[p(z)  j  8o“c*  _  dz _ 

pvt)  “  R  (a  4  *)>  /T  -  U  4  L*\ 

a  o  I 
vo  / 

where  a  -  the  radius  of  the  earth 

gQ  "  the  acceleration  due  to  gravity  at  »ea  level 
=  the  molecular  weight  of  air  at  sea  level 
and  k  =  the  gas  constant 


Then  by  uslrg 


Integration  ot  the  above  between  s  and  *q  remits  in  the  final  form,  which  ia 
.ised  to  interpolate  pressures  in  an  altitude  group 
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To  calculate  the  defining  properties  of  the  model  atmosphere,  the  pressures  at 
the  b#»r?  point  of  each  altitude  group  and  the  molecular  temperature  sea  level 
were  taken  as  constants.  Then,  starting  at  sea  level,  an  iterative  procedure 
was  used  to  f'nd  the  molecular  temperature  gradient,  L  in  the  fir.it  alti¬ 

tude  group.  Using  L,  thus  determined,  the  -cclecular  temperatu* e>  of  the  base 
point  o'  ' next  altitude  group  was  calculated.  Then  the  iterative  procedure 
was  again  used  co  determine  L  for  the  second  altitude  group  With  this  technique 
the  defining  properties  for  the  atmospheric  model  stable  under  the  l/R*  gravity 
field  were  detaraiaed.  The  fincl  results  are  given  in  table  II. 
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Table  11 

DEFINING  PROPERTIES  OF  THE  AFWL  STANDARD  ATMOSPHERE 


Altitude 

z(Km) 

Molecular  seals 
temperature 

T  (*K> 

50 

Gradient 

L  (*K/cm) 

Pressure 
p  (Jyoes/ctn2) 

0.0 

288.1500 

-6.492636 

X 

10“5 

1.01325 

x  106 

11.0 

216.7310 

-1.842590 

X 

10"7 

2.26999 

X  105 

2C.0 

216.5652 

9 .955301 

X 

10"6 

5.52930 

X  104 

32.0 

228,5163 

2-723233 

X 

10"5 

8.89063 

X  103 

47.0 

269,3648 

5.002971 

X 

10"6 

1.15851 

*  103 

52.0 

271.8663 

-2.051281 

X 

1P~5 

6.22283 

X  102 

61.0 

253. 4048 

-3.760658 

X 

10"5 

1.96917 

X  102 

79.0 

135.7129 

-8.818089 

X 

10~6 

1.24420 

-  101 

90.0 

176.0130 

3.936-81 

X 

10"  5 

1.64330 

100,0 

215.8779 

3.882785 

X 

10~5 

3.00750 

X  10"1 

110.0 

254.705”' 

1.134533 

X 

10"4 

7.35440 

X 

*-» 

o 

1 

PO 

120.0 

368.1590 

1.927120 

y 

10"4 

2.52170 

<N 

1 

o 

H 

X 

150  0 

946.2950 

1.803651 

X 

10"4 

5.06170 

X  1(T3 

xuG.O 

1126.6600 

6.716274 

X 

10"H 

3.b9430 

x  10" 3 

170.0 

1193.8230 

8,753717 

X 

10"5 

2 .79260 

x  1G~3 

19C.0 

13b8.8970 

4.U47483 

X 

10"5 

1,68520 

x  10"3 

230.0 

1530.7970 

4.605377 

X 

10"5 

6.96040 

x  10~L 

300.0 

1853.1730 

2.825029 

10“5 

1.88380 

x  10"4 

400.0 

2135.6750 

3.127025 

X 

10"5 

4.03040 

x  10-5 

500.0 

2448.3780 

1.135442 

X 

10”5 

I.C1570 

x  10"5 

600.0 

2561.9230 

1.695913 

X 

10“5 

3.45020 

x  10"6 

700,0 

2731.5140 

•  «  i 

1 

1 .  J.5180 

x  io~6 

28 
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SECTION  V 


RESULTS 


The  output  of  both  SHELL  and  SAP  consists  of  values  for  pressure,  velocity, 
density,  and  specific  Internal  enargy  for  each  zone,  every  cycle.  However,  there 
quantities  are  recorded  on  output  tape  at  selected  tines  only. 

To  present  this  output  in  the  most  compact  form,  the  quantities  on  the 
output  tape  ere  processed  by  a  separate  code  that  generates  plots  of  the  above 
varlabxes.  SAP  plots  consist  of  profiles  only.  SHELL  plots  consist  of  histo¬ 
grams,  contours,  and  velocity  vectors. 

The  overpressures,  dynamic  pressures,  impulses,  and  velocities  recorded  at 
the  "test  stations"  were  also  plotted  as  a  function  of  time  for  each  station. 

Following  normal  procedure,  SHELL  and  SAP  peak  shock  pressures  were  extrapo¬ 
lated  to  find  a  "best"  value  for  these  quantities.  This  is  usually  necessary 
since  numerical  methods  smear  the  shock  over  several  zones.  The  extrapolation 
procedures*  for  SAP  and  SHELL  were  developed  based  on  a  comparison  or  results  of 
SAP  and  SHELL  calculations,  of  a  1-KT  detonation  at  sea  level  with  the  classical 
Taylor  solution  of  a  strong  blast  wave  resulting  from  an  intense  explosion. 
Further  comparisons  in  both  strong  and  weak  shock  regions  were  made  with  the 
1959  blast  curve  (Ref  17),  which  has  been  accepted  to  represent  the  experimental 
data  from  a  1-KT  nuclear  free-air  burst  at  sea  level.  It  was  found  that  SAP  peak 
pressures  should  be  extrapolated  to  that  radius  where  the  maximum  artificial 
viscosity  occurs.  This  procedure  holds  for  both  strong  and  weak  shock  though  for 
strong  shocks  the  extrapolated  pressures  are  not  significantly  different  from 
the  unextrapolated  ones.  It  was  also  found  that  SHELL  peak  pressures  should  be 
extrapolated  to  that  ijdius  corresponding  to  1/2  the  peak  overpressure  for 
strong  shocks  (above  3  atm  overpressure),  but  for  weak  shocks,  the  actual  cal¬ 
culated  peak  pressure  should  be  taken  at  that  radius  corresponding  to  1/2  the 
peak  overpressure.  See  figure  4. 

The  above  procedures  were  used  co  extrapolate  SAP  and  SHELL  peak  overpres¬ 
sures  for  the  methane  and  TNT  detonations.  The  extrapolated  overpressures  are 
plotted  vs.  radius  and  appear  in  figures  5  and  6. 
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Figure  4.  Results  of  Pressure  Extrapolation 
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1.  TNT 

The  following  discussion  of  results  of  the  TNT  calculation  applies  to  a 
scaled  20-short  ton  detonation. 

The  detonation  wave  in  the  SAP-TNT  calculation  reached  a  self-similar  etate 
very  quickly  with  a  peak  pressure  of  1.83  *  1011  dynes/ca2 .  It  reached  the 
surface  of  the  TNT  at  0.2048  Billiseconds.  As  the  shock  encountered  the  air, 
free  expansion  commenced.  The  peak  pressure  began  to  drop  and  a  rarefaction  wave 
started  into  the  detonation  products,  toward  the  center.  Meanwhile,  the  free- 
air  shock  was  building  up.  The  rarefaction  wave  reached  toe  center  at  0.6  Billi¬ 
seconds.  At  this  point,  the  detonation  products  were  "cold"  and  much  more  dense 
than  the  atmosphere  surrounding  thea.  However,  the  detonation  products  had  high 
velocities  and  so  acted  like  «.  piston  pushing  the  atmosphe  e.  A  second  shock 
formed  in  the  detonation  products  just  inside  the  TNT-air  interface.  It  had  a 
velocity  smaller  than  the  velocity  of  the  interface  so  the  second  shock  worked  its 
way  into  the  detonation  products.  However  the  net  motion  ;as  still  outward.  At 
16.3  milliseconds,  the  kinetic  energy  of  the  detonation  pri  ducts  had  dissipated 
to  such  an  extent  that  the  second  shock  began  to  move  inwsii  towards  the  center. 

Meanwhile  the  main  outward  moving  shock  as  calculated  by  SAP,  had  reflected 
off  the  surface  of  the  earth  at  13.36  milliseconds,  with  an  incident  overpressure 
of  140  psi,  and  a  reflected  overpressure  of  730  psi.  SHELL  shows  that  after  the 
main  shock  has  proceeded  along  the  ground  for  approximately  on*  burst  height 
(85  feet),  the  mach  stem  begins  to  form. 

At  34  milliseconds  the  second  shock  reflects  off  the  center  and  proceeds 
outward.  It  reflects  off  the  ground  at  approximately  84  millise'onds  with  a 
reflected  overpressure  of  19  psi.  This  shock  is  very  weak  and  is  quickly  dis¬ 
sipated,  There  are  also  other  very  weak  shocks  which  are  barely  oiscernable. 

The  SHELL-TNT  calculation  indicates  reflection  of  the  main  shock  at  the 
ground  at  13.6  milliseconds  with  a  reflected  peak  overpressure  of  720  pel 
(extrapolated) . 

An  auxiliary  SAP  calculation  was  performed,  using  the  same  SAP  conditions 
which  served  as  input  to  the  SHELL  calculation.  (These  conditions  were  deli¬ 
berately  chosen  at  a  sufficiently  advanced  time  when  the  density  and  internal 
energy  of  the  detonation  product!  were  comoarable  to  that  of  air.)  For  this 
run,  the  equation  of  state  for  air  was  used  in  all  regions:  both  TNT  and  air. 

The  purpose  of  this  auxiliary  calculation  was  to  see  whether  there  would  be 
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significant  dlff trancea,  at  Chest  lsttr  Claes,  between  the  results  of  a  calculation 
the  LSZJC  equation  of  state  for  th*  TNT  zones  and  those  Iron  one  using  the 
•qua t ion  of  state  for  air  for  the  TNT  zones.  The  results  ot  this  calculation 
•bowed  no  discernible  difference  in  peak  pressures  or  arrival  tines.  Therefore, 
going  frow  a  two  equation  of  state  SAP  calculation  to  a  one  equation  of  state 
SHELL  calculation  wee  justified.  See  figure  7. 

2.  Methane 

The  detonation  wave  in  the  SAP-«ethane  calculation  reached  a  self  similar 
state  very  quickly  with  a  peak  pressure  of  3.5  *  107  dynea/co2.  It  reached  the 
surface  of  the  balloon  at  6.3  nilliaeconds. 

As  the  shock  wave  hit  the  surface  of  the  balloon,  a  second  shock  waa  reflected 
toward  the  center  fron  the  relatively  higher  density  sab lent  air.  This  shock 
soon  dissipated  and  the  rarefaction  wave  started  toward  the  center.  The  rare¬ 
faction  wave  reached  the  center  at  18.2  milliseconds.  The  rarefaction  was 
followed  by  a  reflected  shock  which  reached  the  center  at  30  ms. 

Meanwhile  the  main  outward  moving  shock  has  reflected  off  the  ground  at  14.2 
milliseconds  with  an  incident  overpressure  of  7.43  *  106  dynes/ca2  (108  psi) . 

For  this  value  of  Incident  overpressure  theoretical  relationships  give  a  reflected 
overpressure  of  3.77  *  107  dynes/cm2  (550  psi). 

The  second  shock  reached  the  ground  at  62  milliseconds  with  a  reflected  over¬ 
pressure  of  25  psi.  This  shock  is  very  weak  and  quickly  dissipated. 

The  SHELL  methane  calculation  indicated  reflection  of  the  main  shock  at  the 
ground  at  14.4  suLlliseconds  with  a  reflected  peak  overpressure  of  580  psi 
(extrapolated),  'rom  this  time  on,  th°  two-dimensional  phenomenology  is  similar 
to  that  for  the  TNT. 

During  the  methane  calculations  the  peak  overpressures  for  stations  5  and  6 
were  lost.  The  loss  of  these  peaks  affect  only  the  overpressure  waveform  and 
overpressure  impulse  for  these  stations.  All  other  parameters  are  unaffected. 

The  overpressure  plots  which  appear  for  these  stations  (5  and  6)  have  onl those 
points  which  were  dumped  on  output  tape  to  replace  the  missing  points.  The 
appendix  to  this  report  contains  sunsaary  plots  of  information  contained  in  the 
detailed  plots  of  Volume  II. 
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Figure  7.  Comparison  of  Peak  Fraaaure  vs.  Radius  for  the  Equations  of  State  for  TNT  and  Air 
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3.  Comparison  of  the  tjrthm  and  TNT  Results 

Both  runt  wrt  Mule  on  identical  g^lds,  the  only  difference  being  the  input 
conditions  froei  the  SAP  calculation.  In  both  caaea  an  interface  of  trace 
particles  was  placed  at  the  gas-air  Interface. 

Title  interface  brings  to  light  the  first  notable  difference  in  the 
calculation.  The  high  density  TNT  expanding  into  the  lov  density  air  ^ives  rise 
to  Taylor  instabilities  at  the  interface.  These  can  be  clearly  seen  on  the  TNT 
plots.  Since  the  methane-oxygen  mixture  has  nearly  the  sue  density  as  air,  the 
above  condition  does  not  exist,  and  in  fact  is  not  seen  in  the  calculation. 

The  overpressure  for  TNT  ct  a  radius  of  16.8  meters  is  less  than  that  for 
•ethane  at  the  same  radius.  At  a  radius  of  26  meters  the  TNT  overpressure  exceeds 
that  of  the  aethane.  This  is  due  to  the  somewhat  larger  kinetic  energy  of  the 
TNT  shock. 

Although  the  peak  pressure  and  peak  dynamic  pressure  of  the  TNT  shock  are 
greater  than  that  of  methane  to  a  radius  c*  "*4  meters,  the  impulses  in  both 
cases  are  larger  for  the  methane  than  for  TNT.  This  means  that  the  shock  from 
•ethane  is  much  thicker  than  that  of  TNT. 
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SECTION  VI 
CONCLUSIONS 

The  SAP  and  SHELL  codes  have,  In  past  calculations,  reproduced  the  phenome¬ 
nology  of  ataoepheric  blast  both  quantitatively  and  qualitatively.  Results  have 
agreed  with  experimental  data  to  within  10  percent.  Therefore,  there  is  con¬ 
siderable  confidence  in  the  results  reported  here. 

As  can  be  seen  in  figures  8  and  9  there  Is  agreement  between  SAP  and  SHELL 
overpressures,  and  both  are  troa  3  percent  to  20  percent  higher  than  Brode's 
(Ref  16).  The  dynanlc  pressures,  velocities  end  impulses  are  good  to  within 
5  percent.  The  calculated  arrival  time  fer  station  1  agrees  with  Brode's  arrival 
time  to  within  2  oercent. 
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APPENDIX 

This  appendix  presents  a  summary  of  the  calculation  in  graphical  form. 

This  information  was  obtained  from  the  detailed  iistory  of  the  calculation  con¬ 
tained  in  Volume  II. 
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Figure  12.  Ground  Level  Overpressure  Impulse  vs.  Ground  Range. 
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Figure  14.  Dynamic  Pressure  vs.  Ground  Range. 
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